Cracking Piles of Brittle Grains 
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A model which accounts for cracking avalanches in piles of grains subject to external load is 
introduced and numerically simulated. The stress is stochastically transferred from higher layers to 
lower ones. Cracked areas exhibit various morphologies, depending on the degree of randomness in 
the packing and on the ductility of the grains. The external force necessary to continue the cracking 
process is constant in wide range of values of the fraction of already cracked grains. If the grains 
are very brittle, the force fluctuations become periodic in early stages of cracking. Distribution of 
cracking avalanches obeys a power law with exponent r = 2.4 ± 0.1. 

PACS number(s): 64.60.Lx; 81.05.Rm; 46.30.Nz; 83.70.Fn 
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I. INTRODUCTION 

There are many phenomena concerning granular mat- 
ter which attract attention of physicists jlj . The source of 
complexity of sand and similar systems stems from highly 
non-linear mechanic response on the mesoscopic scale (i. 
e. on the scale of single grains) which brings about com- 
plicated behavior on many scales, up to the macroscopic 
one, even though there is usually no scale-free behavior 
H . This feature brings the physics of granular matter 
close to other complex mechanics phenomena, like fric- 
tion H and wear j|, where the interplay of mesoscopic 
and macroscopic phenomena is the central point of at- 
tention. 

The dynamics of sand may be studied from two points 
of view. Slow driving by adding single grains gives rise 
to avalanches and stratification phenomena M. In- 
tense driving by periodic or persistent external forces was 
observed to cause for example surface pattern formation 
(dunes etc.) or grain-size separation jlj. Dynamics of 
mixture of sand and air may lead to beautiful phenom- 
ena like ticking of hour glasses || . 

On the other hand, the most frequently asked question 
about static properties was the stress distribution within 
sand heaps, either free of embedded in various kinds of 
containers pj,p|-|l"3|| . The most famous phenomenon is 
perhaps the minimum of stress directly below the top of a 
conic sandpile, measured by Smid and Novosad jl4] and 
later on explained theoretically by Bouchaud and oth- 
ers p^-|T7||. The explanation is based on the fact, that 
arches are created within the granular packing, which 
support most of the weight. A very important phenom- 
ena connected with arching are the static avalanches due 
to large-scale reconstruction of arches, caused by very 
small external perturbation jl8| , and stick-slip motion of 
sand in a tube fl9| , p0| ]. 

Both of the above phenomena are currently well de- 
scribed within the scalar arching model ]l9p , which is a 
generalization of the scalar stress model developed for 



granular matter by Liu et al. (]l3| , f2l| , [22|| . 

Less studied phenomenon from the point of view of 
granular materials is the procedure in which the grains 
are produced, i. e. the fragmentation process |p3| , p4j| . 
The obvious practical importance of this process was 
stressed e. g. in In the statistical approaches 

to fragmentation p3[ , the grains which are cracked are 
considered either independently of each other or random 
two-particle collisions of the grains are taken into ac- 
count. Such models are appropriate to the situation in 
mills. Different mechanisms should be at work when the 
bulk of the heap of granular particles is cracked by com- 
pression, like in manufacturing pills in pharmaceutical 
industry. Similar problems were already addressed when 
studying localization of deformation in two-dimensional 
heaps of plastic cylinders p^j and compaction of granular 
matter in silos under pressure p7| . 

In the present work, we introduce a model, which con- 
siders cracking of grains within a pile of other grains, 
some of them already cracked, others not. So, we will not 
investigate the size distribution of fragments, like in Ref. 
p3| , but the spatial configuration of clusters of cracked 
grains and also the external force fluctuations occurring 
during the process of cracking. 

The article is organized as follows. In the next section 
the model is introduced. The section III is a gallery of 
simulation results and the last section, sect. IV draws 
conclusions from the results obtained. 



II. DESCRIPTION OF THE MODEL 

Our model describes a two-dimensional pile of granular 
matter contained in a rectangular silo. A physical real- 
ization of this situation may be prepared by two parallel 
glass plates, distance of which corresponds to grain size. 
The lateral and bottom slots are closed, while the upper 
slot is open and a uniform external force is applied to the 
surface of the pile by a kind of piston. The grains are brit- 
tle (eggs may serve as a popular example), which means 
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that if the stress the grain supports exceeds a threshold 
value u>thr, the grain collapses. As a consequence of this, 
the stress pattern in the neighborhood of the collapsed 
grain changes, which may cause another grain collapse 
and finally leads to a kind of internal avalanche. During 
that process, the piston is kept immobile, so the total 
external force decreases, until the avalanche stops. How 
much the force decreases as a consequence of cracking 
one grain, is described by a material dependent factor 
a < 1. We may expect, that for more ductile grains, the 
drop of the force will be smaller and the parameter a will 
be closer to 1. For this reason we will call a the ductility. 

The stress within the pile is a tensor, but recent stud- 
ies [fl7f showed that for many purposes only the diagonal 
element corresponding to the horizontal axis is impor- 
tant. This simplification leads to a scalar model of stress 
propagation in granular matter, which will be a basis of 
our model here. 

We suppose the grains are placed regularly on a square 
lattice rotated by 45 degrees, so that the columns and 
rows of grains correspond to the diagonals on the lattice. 
Each row is L grains wide, each column is H grains high. 
The grains are in contact with the nearest neighbors on 
the lattice. The randomness in the size, shape, and po- 
sition of the grains is taken into account by a stochastic 
rule, which describes the propagation of stress. 

Denote Wik the stress on the grain in z-th row (counted 
from above) and fc-th column. It transfers the fraction 
qik of the stress to its left bottom neighbor, the fraction 
1 — Qik to its right bottom neighbor. We neglect the 
weight of the grains themselves, compared to the external 
force. So, the rule of stress propagation is described by 
the equations 

Wi+l,k = qikWik + (1 - qik-i)wik-i for odd k , . 
Wi+l,k = (1 - qik)mk + Qik+iWik+i for even k . 

We impose cylindrical boundary conditions, wm = WiL- 
The topmost row is subject to external force w\k = fk- 
Total force on the piston is then F = J^k fk- 

The simulation proceeds as follows. The numbers qik 
are taken randomly from the uniform distribution on the 
interval (J-^-, ^~^)- Initially all fk are set equal and the 
local stresses are computed according to rules (|l|). The 
force is increased until stress on one non-cracked grain, 
say at position (i, k), reaches the threshold w t h r = 1. 
Then, the grain is cracked, which has two consequences. 
First, the force on top of its column is lowered, fk — > afu- 
Then, if grain in the same row to the left, i e. (i, k — 1) 
is not cracked, the value of q corresponding to left top 
neighbor of (i, k) is set to 1. If (i, k — 1) is cracked q is 
given new random value from the uniform distribution 
on the interval (^^, ^~2^)- Similar rule applies on the 
right hand side: if (i, k + 1) is not cracked, the right top 
neighbor if (i, k) has new q = 0, if (i, k + 1) is cracked, 
the new q is a random number from the same distribution 



as above. These rules correspond to very simple intuitive 
observation, that the cracked grain does no more bear the 
load, if it has neighbors, which can bear the load instead 
of it. However, if the neighbors are also already cracked 
grains, the stress propagation remains to be stochastic 
as it was before the cracking, but the realization of the 
randomness, i. e. the values of the numbers q are 
changed. 

After each change of q's, the local stresses are recom- 
puted, the grains which are not yet cracked and exceed 
the threshold are cracked, new q's are established and 
this procedure is repeated until no non-cracked grains 
exceeding the threshold are found. Then, the external 
force is increased up to the value when another grain is 
cracked again and new cracking avalanche begins. We 
will call avalanche size s total number of grains cracked 
during the avalanche. This algorithm continues as long 
as there are any non-cracked grains left. 

Besides the size of the system, the model has two free 
parameters. The parameter a measures the ductility of 
the grains and (3 the degree of randomness in the stress 
propagation. The limit case j3 = corresponds to fully 
deterministic case. 

III. SIMULATION RESULTS 

When a grain is cracked, the load is mostly transferred 
to its neighbors, which have then increased chance to be 
cracked. This leads to creation of clusters of cracked 
grains, which grow and merge as the cracking proceeds. 
Typical morphology of the cracked clusters is shown in 
the Fig. 0. We can observe formation of "arches" with 
one dominant "leg" only. The shape of the "legs" resem- 
bles the letter S when they grow large. The dependence 
of the morphology on the ductility a and randomness /3 
is shown in the Figs. 0, ||, and ||. For larger (3 the typi- 
cal size of the cracked clusters is smaller, while for small 
(3 the sample contains only few big "arches" , which are 
also more symmetric than those for larger randomness. 
The ductility has different influence on the morphology: 
in the case of more brittle grains, i. e. with smaller a, 
the cracked areas are mostly concentrated in the top part 
of the sample, while more ductile grains lead to cracking 
equally probable in the whole bulk of the sample. (We 
performed simulations also for very ductile grains, a close 
to 1, and the trend was observed to shift the cracked re- 
gions to the bottom of the sample, when the ductility is 
increased.) 
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FIG. 1. Morphology of cracked areas for a sample with 
L = 500, H = 500, after 5000 time steps. Every cracked grain 
is represented by a black dot. The parameters are a = 0.9, 
(3 = 0.25. 
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FIG. 3. Morphology of cracked areas for a sample with 
L = 500, H — 500, after 5000 time steps. The parameters are 
a = 0.9, /3 = 0.1. 
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FIG. 2. Morphology of cracked areas for a sample with 
L = 500, H = 500, after 5000 time steps. The parameters are 
a = 0.9, (3 = 0.5. 



FIG. 4. Morphology of cracked areas for a sample with 
L = 500, H = 500, after 10000 time steps. The parameters 
are a = 0.5, (3 = 0.5. 
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When the cracking proceeds, the force necessary to 
continue fluctuates. Each cracking avalanche means a 
drop of the force, which then rises again. The Fig. | 
shows the time dependence of the external force F and 
the fraction of cracked grains v for the sample of 200 x 200 
grains. We can see that the force fluctuates around a 
nearly time-independent value F w ~ 0.55 during large 
part of the process, at least from time t = 1000 to 
t = 5000. This was even more clearly observed for larger 
samples (in our simulations 500 x 500). So, the picture of 
the overall behavior of the force can be as follows. After 
a transient period, where the force suddenly drops and 
slowly rises again, a stationary cracking regime devel- 
ops, characterized by constant average force F av . This 
regime holds if the fraction of cracked grains is small, 
according to our observations v < u max is sufficient con- 
dition, where the value of z/ max depends slightly on a. 
For a = 0.1 we found ^ max — 0.7, while for a — 0.9 we 
observed z^ max ~ 0.4. 



with a power-law tail seems to be a good candidate, but 
further data would be needed to settle this question. 
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FIG. 5. Time evolution of external force F (full line) and 
fraction of cracked grains v (dashed line) for the sample with 
L = 200, H = 200, a = 0.9, /3 = 0.25. 

The value of stationary force F av decreases with j3. 
we found the values in the range from F av ~ 0.3 for 
(3 = 1 (maximum randomness) to F av ~ 0.6 for [3 = 0.1 
(minimum randomness studied). 

Around the average force, there are fluctuations, which 
reflect unique realization of the disorder in our sample. 
We investigated statistical properties of the fluctuations, 
using histogram of the changes of force from one time step 
to the next one. The distribution of upward changes can 
be very well fitted by an exponential, while the down- 
ward changes do not have any clear form of distribu- 
tion: neither Gaussian, exponential, stretched exponen- 
tial nor power-law fit was satisfactory. A distribution 
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FIG. 6. Time evolution of external force F (full line) and 
fraction of cracked grains v (dashed line) for the sample with 
L = 200, H = 200, ol = 0.1, (3 = 0.25. 

For very small a (we observed the phenomenon for 
a = 0.1, but for a — 0.3 it was already absent) the 
fluctuations loose their purely random appearance and 
quasi-regular force oscillations occur, which are especially 
pronounced in the early stages of the cracking process (i. 
e. for small v). They can be clearly seen in the Fig. 
||. When the fraction of cracked grains increases, the 
oscillations gradually disappear. The oscillations perhaps 
correspond to the sudden drop of the force, observed for 
all a, followed by gradual increase of the force again. 
While for small a many periods of the oscillation may be 
realized, for larger a the oscillations are "over-damped" 
and only single period occurs. 

A cracking avalanche starts from the stable state, in 
which stresses on all non-cracked grains is below the 
threshold. The avalanche is initiated by increase of ex- 
ternal force up to value which causes one grain to crack. 
This cracking may result in cracking other grains, and so 
on, until new stable state is reached and the avalanche 
stops. We denote Ac the avalanche size, which is the 
number of grains cracked during the avalanche. We are 
interested in statistical distribution of avalanche sizes. 
We expect, that the distribution may be different in the 
initial transient period and in the stationary regime, in 
which the average force f av is constant. So, we investi- 
gated the distributions P i > (Ac) defined as probabilities 
that the size of the avalanche, occurring in time interval 
(t/_i,ij], with to = 0, is larger than Ac. 
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Figure [?] shows the results for a 500 x 500 sample. The 
first two intervals, with final times t\ — 100 and = 300 
describe the situation in the transition regime. We can 
see that most of avalanches have typical size about Ac ~ 
400. On the other hand, the next two intervals with 
end times £3 = 1000 and £4 = 5000 give distributions 
which can be fitted by a power law in the range of two 
decades. It can be also seen that the distribution is stable 
in time during the stationary cracking regime. We fitted 
the exponent of the power-law dependence P > (Ac) ~ 
(Ac) 1_r with the result 

t= 2.4 ±0.1 . (2) 

We have found the same exponent (within error bars) 
for all values of parameters studied. The only exception 
was the case of j3 = 1, where the distribution was close 
to exponential, instead of power-law. The breakdown of 
power-law, when (3 approaches 1 remains to be studied. 




FIG. 7. Avalanche size distributions for L — 500, H = 500, 
a = 0.9, (3 = 0.25, in intervals determined by times t\ — 100, 
ti = 300, t 3 = 1000, and U = 5000. The lines denote the 
following distributions: dash-dotted line P{* , dotted line P 2 > 1 
dashed line P^" , solid line P± . P{* corresponds to the interval 
(ti-i,ti\. 

IV. CONCLUSIONS 

We have found that the two-dimensional pile of brit- 
tle grains packed in a rectangular container exhibit non- 
trivial behavior, when an external force is applied from 
above and the grains are cracked. The cracked grains 
form clusters with different morphologies, depending on 
the ductility of the grains and on the degree of random- 
ness in the packing. Degree of randomness seems only to 
determine the characteristic scale of the cracked clusters: 



lower randomness leads to larger clusters. This fact can 
be understand rather easily, if we realize that a cluster 
occurs when the local stress exceeds the threshold neces- 
sary for a grain to be cracked. If the stress distribution is 
more uniform, fluctuations above the threshold are more 
distant one from the other. 

Less expected feature is the influence of the ductility. 
Brittle grains have tendency to crack in the top part of 
the container, while ductile grains are cracked mostly in 
the bottom part. This finding may play important role 
in separation of grains of different types. 

During the cracking process the external force fluctu- 
ates around a general trend, which can be described as 
follows. If the grains are not too brittle (a > 0.3), the 
force drops suddenly and then rises slowly to a value, 
which then remains constant for great part of the whole 
cracking process. When the fraction of cracked grains 
approaches 1, the force increases again. So, there exists 
well-defined stationary cracking regime, preceded by a 
transient period and followed by a final stage. For very 
brittle grains, the force oscillates rather regularly even in 
the stationary regime. 

Cracking one grain may result in an avalanche of fur- 
ther crackings. The distribution of avalanche sizes de- 
pends on time. While in the transient period the distri- 
bution is not scale-invariant, in the stationary regime the 
distribution of avalanche sizes obeys a power law. This 
is an indication, that a sort of criticality is present in 
the cracking process. The value of the exponent r ~ 2.4 
is larger than avalanche exponents found in most self- 
organized critical models known to us. On the other 
hand, the dynamics of our model resembles the Olami- 
Feder-Christensen model of earthquakes p8| , where the 
exponent varies in wide range, comprising also the value 
found in our model. However, the mechanism leading 
to power-law scaling in the OFC model is not completely 
clear. This may suggest that a new mechanism leading to 
criticality is at work here, different from the usual SOC. 

This work does not compare the simulation results with 
experimental data, because we were not able to find any 
report of an experiment of this kind (loosely related are 
the experiments reported in fffi] ). It would be very wel- 
come if a measurement in the direction suggested here 
was done in future. 
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